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Abstract. We investigate bifurcation phenomena between slow and fast convergences 
of synchronization errors arising in the proposed synchronization system consisting of 
two identical nonlinear dynamical systems linked by a common noisy input only. The 
numerical continuation of the saddle-node bifurcation set of the primary resonance of 
moments provides an effective identifier of the slow convergence of synchronization er- 
rors. 
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1. Introduction. The noise-induced synchronization of a dynamical system with its 
copies can easily be found in nonlinear systems, such as the discrete mapsP, [2], the 
Lorenz system[T], the Duffing oscillator [3], the single mode CO2 laser and the un- 
coupled neurons[S]. One of the most important results of them is that the perfect syn- 
chronization may arise under some suitable conditions^, El II]. Moreover, the perfect 
synchronization exhibits significant degree of robustness against mismatches among the 
copies such as the parameters mismatch [3] and the independent random fiuctuations of 
the copies [H [6]. Furthermore, regarding the response of the synchronization system as 
a Markov process to derive the transition law of it, we have analytically shown that the 
perfect synchronization can be regarded as an absorbing barrier of the Markov process[7]. 

In these studies, however, little research have been done on transient behavior to ap- 
proach the perfect synchronization. In engineering applications, too slow convergence of 
synchronization errors would be regarded as failing to converge, even if the synchronization 
is achieved mathematically. In practice, the perfect synchronization is possibly applicable 
to synchronizing initial conditions of coexisting oscillators with a common specification 
such as independent subcircuits in a circuit system, independent mechanical vibrators on 
bench testing, and so on. However, it is hardly applicable to industrial purposes when 
the convergence speed is too slow. 

To solve this problem, we have already investigated how to characterize the slow con- 
vergence of synchronization errorsjHl E], showing that the slow convergence is related to 
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the slow diffusions caused by multimodal probability densities so that it can be detected 
as multi- valued solutions of the moment differential equations (MDE). 

In this paper, we perform nonlinear analysis on the MDE. We first construct Poincare 
maps of periodic solutions of the MDE to examine asymptotic behavior of the moments 
and clarify that the primary resonance encounters a saddle-node bifurcation. We then 
numerically continue the bifurcation point to obtain the saddle-node bifurcation set that 
is in good agreement with the threshold of occurrence of the slow convergence. 

2. Noise-induced Synchronization. 

2.1. Slow convergence of synchronization errors. Let us consider the pair of iden- 
tical piecewise linear systems[8j, |9j: 

X + cx + kG{x; ji) = Q + u{t), , . 

y + cy + kG{y; fi) =Q + u{t). ^ ^ 

where Q is a preload, u{t) is a random input, and the function G is a piecewise linear 
function defined hj G{x; fi) = {x + fi) — \x + fi\ + {x — fj,) + \x — fi\ describing a linear 
spring with a dead zone of the width 2/i. In the previous work[D], we have examined the 
two types of the input u{t) such as the combination of the harmonic forcing and white 
noise: 

u{t) = P cos(cut) + sw{t), (2) 

and the filtered noise: 

U + 2C,UJnU + UJnU = Sw{t), (3) 

where w{t) is the standard Gaussian white noise. It is shown that the synchronization 
system ([T]) produces the perfect synchronization where the synchronization error vanishes 
deterministically. 

On the contrary, in the present work, we consider the random-phase harmonic forcing 
as a purely random input in the following form: 

u{t) = P cos [ujt + pBi^ (4) 

where Bt is the standard Brownian motion. We choose c = 0.04, k = 1.0, fi = 0.7, 
Q = 0.3, P = 0.2, and p = 2 x 10~^ which produce both the fast and slow convergence 
by sweeping uj and can be implemented easily as a realistic mechanical structure. 

Figure [T] shows sample paths of synchronization errors for u = 0.8 (upper side) and 
0.95 (lower side) respectively. There is significant difference in convergence time between 
the two conditions. The convergence times are t > 20000 and t = 118 respectively. This 
means that one spends time to converge more than 169 times as much as the other. 

2.2. Localization of probability densities. One explanation of the slow convergence 
of the synchronization errors is the transient localization of probability densities as already 
discussed by the authors [9]. 

One of the pair ([T]), say a single system, produces the transient probability density 
p{xi,X2,t) satisfying the following Fokker-Plank-Kolmogorov (FPK) equation: 



dp d^p d{ujp) d{x2p) d 
dt 2 dxl 8x3 dxi dx^ 



— CX2 — kG{xi; fi) + Q + P cosxs^p (5) 
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Figure 1. Sample paths of synchronization errors for uj = 0.8 (upper) and 
0.95 (fower) respectively. 



where Xi := x, X2 ■= x, X3 := Bf. To obtain direct solutions of these equations, we employ 
Monte-Carlo techniques where the densities are obtained from cumulative frequencies over 
the 10000 samples of numerical solutions of one of the pair ([T]). 

Figure [2] shows dependency of the solution p{xi,X2,t) of the FPK equation ([5]) on 
the initial conditions p{xi,X2,0). The small graphs represent the probability densities 
at t = 60 X 2tt/uj starting from the initial functions p{xi,X2,0) = 6{xi — 2)6{x2 — 2) in 
the upper side and 6{xi + 0.5)5(x2) in the lower side respectively. It is obviously seen 
that the transient densities significantly depend on the initial conditions within the region 
UJ G [0.76,0.92]. In case of u; = 0.76, both the initial conditions (a) and (b) cause quite 
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Figure 2. Transient isolation of probability densities at t = 60 x 2tt/u 
comparable with deterministic hysteric jumps. 



similar densities so that dependency on the initial conditions is hardly found in this case. 
As u increases to 0.78, the outer ring found at a; = 0.76 is being replaced with the smaller 
ring for (a) while the previous state maintains for (b). Further increase of u cause second 
change of the densities at = 0.92 where the outer ring of (b) ai uj = 0.9 is being 
replaced with the smaller ring of (b) at u = 0.92 while the densities of (a) maintains 
between uj = 0.9 and 0.92. 

In the previous work [9], we have already discussed that this dependency on initial 
conditions, in other words, the transient isolation of the densities is one of the main 
reasons of the slow convergence of the synchronization errors and can be detected by 
hysteric jumps of statistical moments. It will be shown that similar results can be obtained 
even in this work considering the random-phase forcing (jl]). 

2.3. Moment equations. In order to evaluate the statistical moments of the synchro- 
nization, we derive the moment differential equations (MDE) from the FPK equation ([5]) 
as follows. 

nil = fn2, = —cm2 — k{G) + Q + P(cosx3), 
wis = tu, a'li = 2(Ti2, 

0"12 = 0"22 - C(Ti2 - k{{xi - mi)G) + P{Xi COSX3), 

Cr'l3 = 0^23 + i^mi, 

(722 = -2C(T22 - 2A;((X2 - 1712)0) + 2P(x2COSX3), 
0-23 = -CO-23 - k{{x3 - ma)^) + P{X3 COSX3) - 001712, 

(T33 = 2ujm3 + (6) 

where mj = (xj), 0"^ = ((xj — mi){xj — TUj)). The nonlinear averages are linearized by 
the standard statistically equivalent techniques pUl [TT] . In our case, however, 2:3 = Bt is 
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hardly assumed to be Gaussian because is unbounded. Therefore, utihzing the formula 
of the Brownian motion [I^: 

1 f°° 

U{Bt)) = ^l f{x)e-^'/'dx, (7) 



oo 



we perform the approximation (cosxs) ~ /^cos(m3), (0:3 cos X3) ~ ^'^ ^/ujrn^ cos m^. 

In this setup, the moment equations © can be regarded as a periodically forced system 
(reducing fn^ = a; to = ujt) where the third equation of ^ is neglectable. Furthermore, 
in order to calculate an and (T22 only, the sixth, eighth, and ninth equations are also 
neglectable because the remaining equations are not coupled with them. Therefore, we 
can reduce the MDE ([6]) into the following form: 

nil = "^2, 

7712 = — cm2 — kao{mi, an) + Q + Pe^ cos cut, 
a'li = 2(Ti2, 

0"i2 = cr22 — cai2 — kai{mi,aii)aii + Prriie''^^'^ cosut, 
(J22 = — 2co"22 — 2kai{mi, an)o'i2 + 2Pm2e^^^^ cos cut (8) 
where ao{mi, au) and ai(mi, an) are the statistically equivalent gains derived by Sato[llJ: 



ao = mi + xl—\e +e ^^rn 

a, = 1 + lerf _ ^ rf f-!^") (10) 



where erf(-) is the error function. Therefore, in what follows, we calculate an and (J22 
from final expression ([8]). 

2.4. Jumps of moments. Figure [3] shows the peak-to-peak values of {xl) = a22 calcu- 
lated from the reduced MDE ([8]) and the mean convergence time [9]: 

(T):=lf;min|n sup |X^(n) - X^(n)| < e 1 (11) 

where n is a discrete time s.t. t = nAt, {X™(n)}^^;^ is the fc-th set of M points com- 
posing numerical probability densities, k represents the fc-th sample path of the random 
excitation, and e <^ 1 is a criterion of the convergence. The moments are calculated as 
numerical solutions of the moment equations ([H]). The solid and broken curves indicate 
the forward and backward sweeps of the frequency u respectively. The mean convergence 
time (T) is estimated by substituting numerical solutions of ([1]) starting from 5x5 uni- 
form Cartesian grids on the region {xi,X2) G [—10,10] x [—10,10] into Eq.l ITT]) where 
M = 5 X 5 = 25, K = 100, e = 10~^. The plots of (T) are saturated at the given value 
(T)„iax = 3000. 

It is clearly seen from Fig 13] that the mean convergence time (T) rapidly increase within 
the hysteric jumps of the second moment of response {xf). This means that the hysteric 
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Figure 3. The peak-to-peak of the second moment (x^)p_p and the mean 
convergence time (T) of the pair ([1]). 

jumps of the moment in the statistical equivalent sense act as an identifier to indicate the 
slow convergence of the synchronization errors. This means that the present result based 
on the random phase forcing in Eq. (jl]) agrees with the previous result |^ considering the 
other types of forcing in Eqs ([2]) and (l3l). 

3. Bifurcation Analysis of Moments. 

3.1. Poincare map of moments. Let us now rewrite the MDE ([8]) in the vector form: 

q = fiq,t;p), qiO) = qo (12) 

where t E R is time, p G i? is a free parameter, and q = (mi, m2, an, cri2, 0-22) '■ R R^, 
f : R^ X R ^ R^. The solution of this initial value problem can be written in the form: 



(13) 



where is the shift operator generated by the MDE. In our case, the function / is 
supposed to be periodic in time, 

f{q,t + T;p) = f{q,t;p) (14) 

where r = 27r/ti; is a period of the external forcing of the MDE (fT2|) . Now, we provide a 
brief summary of how to characterize periodic solutions of the MDE. See the reference 
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Figure 4. Poincare plots of solutions of the MDE for to = 0.9 starting 
from q{0) = a= (0.5, 2.5, 10, 2, 8) and 0.98a. 

for details. In order to reduce dimension, the Poincare map is defined by. 

In practice, a series of points generated by the Poincare map is obtained by q^. = T'^^qo) = 
(p{kT, ^q) (/c = 1, 2, ■ • • ). Plotting qr^ is sometimes referred to as the Poincare plot. 
Therefore, the m-periodic solutions passing through q E satisfies 

T^{q) = q, (16) 

so that the periodic solutions can be regarded as the fixed points of the mapping T. 

Figure H] shows Poincare plots of the solutions of the MDE ([8]) for u = 0.9, which 
is inside of the hysteric region of the response in Figj3l starting from the initial points 
q{0) = a = (0.5, 2.5, 10, 2, 8) and q{0) = 0.98a. It is clearly seen that the initial point 
q{0) = a converges the 1-periodic attractor far from the origin while q{0) = 0.98a 
converges to a 29-periodic attractor near the origin. In this case, it seems that there is 
a repeller near the initial points so that slight difference in initial conditions results in 
distinct amplitudes of response as already shown in FiglSl 

On the other hand, FigE] shows Poincare plots for u = 0.95, which is outside of the 
hysteric region of the response in FigOl starting from the same initial conditions as those 
in FigJH It is obviously shown that the attractor far from the origin vanishes and the 
periodic attractor near the origin is replaced with a quasi-periodic attractor corresponding 
to a solution with infinite period. Figure [6] shows that stationary probability density of the 
original equation ([T]) under the condition corresponding to the quasi-periodic attractor in 
Figl5]for uj = 0.95. The two peaks of the density are caused by temporal switching orbits 
of the sample paths fluctuating around bistable periodic orbits of the original equation 
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Figure 5. Poincare plots of solutions of the MDE ([8]) for = 0.95 starting 
from q(0) = a= (0.5, 2.5, 10, 2, 8) and 0.98a. 




Figure 6. Stationary probability density of the original SDE ([T]) under 
the condition corresponding to the quasi-periodic attractor in Figj5]for uj = 
0.95. Quasi-periodicity is hardly observed. 



([T]) under the deterministic limit. It follows that the quasi-periodicity of the MDE does 
not reflect this bistability of the sample paths. One explanation of this disagreement is 
the approximation errors of the nonlinear moments. 
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Figure 7. Characteristic polynomial P{s) = \sl — DT{q)\ for s = 1 as a 
function of a free parameter uj. -P(l) vanishes at the saddle-node bifurcation 
point. 

From the above result, we can suppose that the primary resonance (the large amplitude) 
in the response curve in FigjS] undergoes a saddle-node bifurcation at a; ~ 0.91 because 
the pair of the repeller and the 1-periodic attractor in FiglUfor u = 0.9 vanishes in Figl5] 
for u = 0.95 while the off resonance undergoes other types of bifurcation such as pitchfork 
bifurcations, Hopf bifurcations, and so on. 

In what follows, we focus on the primary resonance undergoing the saddle-node bifur- 
cation near u = 0.91 to investigate how it indicates the slow convergence of the synchro- 
nization errors. 

3.2. Numerical continuation of bifurcation points. Let us summarize how to obtain 
bifurcation points numerically. See Refs[T3], [Hj for more details. 

It is mathematically proved that a matrix representation DT{q) of the linearization of 
the Poincare map T of the ordinary differential equation ffT^ is given by, 

DT{q) := J{r), j = Df{q, t; A) J, J(0) = / (17) 

where Df{q,t; A) is the Jacobian matrix and I is the unit matrix. It can also be shown 
that the stabilities of the fixed point q of the Poincare map T is characterized by the 
characteristic value s satisfying the characteristic equation: 

P{s) = \sl - DT{q)\ = 0. (18) 

From the basic theory of bifurcation [ll], it is known that saddle-node bifurcations occurs 
at s = 1. 

Figure [7] represents the value of characteristic polynomial for s = 1 calculated from Eqs 
(fT^ and (|T7|1 for c = 0.04. It is clearly seen that the value -P(l) vanishes near u = 0.91. 
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Figure 8. Numerical continuation of the saddle- node bifurcation point of 
periodic solutions of ([H]) and mean convergence time (T) of the synchro- 
nization errors. The white circles are plotted for (T) < 10'^ and the black 
circles for otherwise respectively. 



Since zeros of characteristic polynomials generally indicate bifurcation points |jl3j, it is 
numerically proved that the primary resonance of the second moment undergoes the 
saddle-node bifurcation near uj = 0.91. 

3.3. Saddle-node bifurcation sets indicating the slow convergence. In order to 
continue the bifurcation points in the parameter space [H], we define new state vectors 
p .= {q, \) G and solve the equations (ITUI) and (|TS|) simultaneously: 

T"^iq) = q, Pis) = \sI-DT{q)\=0. (19) 

The solution p is called a bifurcation set of the map ffTUj) . and consequently yields a 
bifurcation set of the corresponding periodic solution of the MDE (|H]). 

Figure [H] shows the saddle- node bifurcation set in the (a;,c)-plane of periodic solutions 
of the MDE ([8]) and mean convergence time (T) of the synchronization errors. In practice, 
we consider the saddle-node bifurcation point in FigIT] for c = 0.04 as the initial point 
and regard u as the bifurcation parameter A. We then solve the simultaneous equations 
( |T9l) by standard Newton's methods. Shifting the second parameter as c = 0.04 + j ■ Ac 
(j = 1, 2, ■ ■ ■ ), we obtain a series of p on the (a;, c)-plane where Ac is a small value chosen 
for a good convergence of the Newton's methods. 
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The solid curve indicates the saddle- node bifurcation set of the MDE (|H]). The white 
and the black circles represent the fast (T) < 10^ and the slow (T) > 10^ convergences 
where (T) is calculated from Monte-Carlo simulations of the synchronization system (Q. 
Along the bifurcation set, the solid curve, the characteristic polynomial takes zero values 
at which the stability changes. In other words, the changes of parameters crossing this 
curve transversely results in the stability change of the system. 

It is very clear from the above result that the saddle-node bifurcation set belonging to 
the primary resonance indicated by the solid curve in FiglH] is in good agreement with the 
threshold of occurrence of the slow convergence. 

4. Conclusions. We have investigated how to characterize the slow convergence of the 
synchronization errors of the synchronization system which consists of the pair of the 
piecewise linear systems subjected to the common random excitation. We first have 
demonstrated numerically that the convergence speed of the synchronization errors sig- 
nificantly depends on the parameter conditions. It is shown that the slow convergence is 
caused by the multimodal probability densities so that it can be detected as the multi- 
valued solutions of the moment differential equations. 

We then have constructed Poincare maps of periodic solutions of the MDE to examine 
asymptotic behavior of the moments and clarify that the primary resonance encounters 
saddle-node bifurcations. We have numerically continued the bifurcation point to obtain 
the saddle-node bifurcation set of primary resonance of moments. It is clearly shown that 
the saddle-node bifurcation set belonging to the primary resonance is in good agreement 
with the threshold of occurrence of the slow convergence. 

The above result leads to the conclusion that the saddle-node bifurcation set of the 
moment equations provides an effective identifier to detect the slow convergence of the 
synchronization errors. 
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